#WIC co-op eWIC new diff in diff pictures
#last modified: 6 january 2025
#last modified by: Charlotte Ambrozek
#this do file takes estimates from the csdid2 packagefor groups, years, and subsamples and constructs graphs as needed
library(haven)
library(tidyverse)
library(RColorBrewer)
library(ggpmisc)

event <- read_dta("./analysis/output/csdid_tip_event.dta") 

picture <- filter(event, state == "all") %>%
  group_by(outcome, subsample, controls, ev_yr) %>%
  mutate(subsample = gsub("all","All", subsample)) %>%
  mutate(subsample = gsub("chain","Chain", subsample)) %>%
  mutate(subsample = gsub("indep","Independent", subsample)) %>%
  mutate(subsample = gsub("highwicelig","High WIC eligibility", subsample)) %>%
  mutate(controls = gsub("no", "No", controls)) %>%
  mutate(controls = gsub("yes", "Yes", controls))

# make three pictures:
#one that just shows all
#one that compares chain and indep
#one that compares all to highwicelig

setwd("./analysis/output/graphs")

all <- filter(picture, subsample == "All" & controls == "No")
allpre <- filter(picture, subsample == "All" & controls == "No", ev_yr < 0)
allpost <- filter(picture, subsample == "All" & controls == "No", ev_yr >= 0)

all_plot <- ggplot(all, aes(x = ev_yr, y = b)) +
    geom_linerange(aes(ymin = ll, ymax = ul),
                   alpha = 0.7,
                   size = 0.8) +
    geom_point(size = 3) +
    geom_vline(xintercept = 0, lty = 3, colour = "red") + 
    geom_hline(yintercept = 0, lty = 2) +
    geom_line(data = allpre, aes(x = ev_yr, y = avgb), colour = "gray70") + 
    geom_ribbon(data = allpre, aes(ymin = avgll, ymax = avgul),
                   alpha = 0.25,
                   colour = "gray90") +
    geom_line(data = allpost, aes(x = ev_yr, y = avgb), colour = "gray70") + 
    geom_ribbon(data = allpost, aes(ymin = avgll, ymax = avgul),
                   alpha = 0.25,
                   colour = "gray90") +
    labs(y = "ATT estimate",
    x = "Event Year",
    caption = paste0("ATT estimate for the post period shown here is ",
                     as.character(round(allpost$avgb[1], digits = 3)),
                     " with a 95% confidence interval of\n(",
                     as.character(round(allpost$avgll[1], digits = 3)),
                     ", ",
                     as.character(round(allpost$avgul[1], digits = 3)),
                     ").")) +
    theme_light() +
    scale_x_continuous(breaks = seq(-4, 4, 2)) +
    theme(plot.caption = element_text(hjust = 0))
  
ggsave("csdid_tip_auth_all.png", plot = all_plot, width = 6.5, height = 4, units = "in")

allcontrols <- filter(picture, subsample == "All" & controls == "Yes")
allcontrolspre <- filter(picture, subsample == "All" & controls == "Yes", ev_yr < 0)
allcontrolspost <- filter(picture, subsample == "All" & controls == "Yes", ev_yr >= 0)

allcontrols_plot <- ggplot(allcontrols, aes(x = ev_yr, y = b)) +
    geom_linerange(aes(ymin = ll, ymax = ul),
                   alpha = 0.7,
                   size = 0.8) +
    geom_point(size = 3) +
    geom_vline(xintercept = 0, lty = 3, colour = "red") + 
    geom_hline(yintercept = 0, lty = 2) +
    geom_line(data = allcontrolspre, aes(x = ev_yr, y = avgb), colour = "gray70") + 
    geom_ribbon(data = allcontrolspre, aes(ymin = avgll, ymax = avgul),
                   alpha = 0.25,
                   colour = "gray90") +
    geom_line(data = allcontrolspost, aes(x = ev_yr, y = avgb), colour = "gray70") + 
    geom_ribbon(data = allcontrolspost, aes(ymin = avgll, ymax = avgul),
                   alpha = 0.25,
                   colour = "gray90") +
    labs(y = "ATT estimate",
    x = "Event Year",
    caption = paste0("ATT estimate for the post period shown here is ",
                     as.character(round(allcontrolspost$avgb[1], digits = 3)),
                     " with a 95% confidence interval of\n(",
                     as.character(round(allcontrolspost$avgll[1], digits = 3)),
                     ", ",
                     as.character(round(allcontrolspost$avgul[1], digits = 3)),
                     ").\nCounty level controls - median income, population, population < 5, share Hispanic,\nshare Black, share earning < FPL, share earning <180% FPL, and share receiving\nSNAP or cash benefits.")) +
    theme_light() +
    scale_x_continuous(breaks = seq(-4, 4, 2)) +
    theme(plot.caption = element_text(hjust = 0))
  
ggsave("csdid_tip_auth_all_controls.png", plot = allcontrols_plot, width = 6.5, height = 4, units = "in")

# Define the colors for each chain status
colors <- brewer.pal(3, "Set2")[1:2]

# Define the shapes for each chain status
shapes <- c(16, 17)

chainindep <- filter(picture, (subsample == "Chain" | subsample == "Independent") & controls == "No") %>%
    group_by(outcome, subsample, ev_yr) 
# Add new columns for color and shape mapping
chainindep$color <- factor(chainindep$subsample)
chainindep$shape <- factor(chainindep$subsample)
chainindeppre <- filter(chainindep, ev_yr < 0) %>%
    group_by(outcome, subsample, ev_yr) 
chainindeppost <- filter(chainindep, ev_yr >= 0) %>%
    group_by(outcome, subsample, ev_yr) 

# Create the chain/independent plot
chainindep_plot <- ggplot(data = chainindep, aes(x = ev_yr, y = b, color = color, shape = shape)) +
    geom_linerange(aes(ymin = ll, ymax = ul),
                   alpha = 0.7,
                   size = 0.8) +
    geom_point(size = 3) +
    geom_vline(xintercept = 0, lty = 3, colour = "red") + 
    geom_hline(yintercept = 0, lty = 2) +
    scale_color_manual(values = colors, name = "Chain status") +
    scale_shape_manual(values = shapes, name = "Chain status") +
    geom_line(data = chainindeppre, aes(y = avgb, x = ev_yr, color = color)) + 
    geom_ribbon(data = chainindeppre, aes(ymin = avgll, ymax = avgul, color = color),
                   alpha = 0.1) +
    geom_line(data = chainindeppost, aes(y = avgb, x = ev_yr, color = color)) + 
    geom_ribbon(data = chainindeppost, aes(ymin = avgll, ymax = avgul, color = color),
                   alpha = 0.1) +
    labs(y = "ATT estimate", x = "Event Year", 
    caption = paste0("ATT estimate for the chain subsample in the post period is ",
                as.character(round(filter(chainindeppost, subsample == "Chain")$avgb[1] , digits = 3)),
                " with a 95% confidence\ninterval of (",
                as.character(round(filter(chainindeppost, subsample == "Chain")$avgll[1], digits = 3)),
                ", ",
                as.character(round(filter(chainindeppost, subsample == "Chain")$avgul[1], digits = 3)),
                ").\nATT estimate for the independent subsample in the post period is ",
                as.character(round(filter(chainindeppost, subsample == "Independent")$avgb[1], digits = 3)),
                " with a 95% confidence\ninterval of (",
                as.character(round(filter(chainindeppost, subsample == "Independent")$avgll[1], digits = 3)),
                ", ",
                as.character(round(filter(chainindeppost, subsample == "Independent")$avgul[1], digits = 3)),
                ").")) +
    theme_light() +
    scale_x_continuous(breaks = seq(-4, 4, 2)) +
    theme(plot.caption = element_text(hjust = 0), legend.position="none") + 
    facet_wrap(~ subsample, ncol = 2)
  
ggsave("csdid_tip_auth_chainindep.png", plot = chainindep_plot, width = 6.5, height = 4, units = "in")

# make chain independent plot with controls
chainindepcontrols <- filter(picture, (subsample == "Chain" | subsample == "Independent") & controls == "Yes") %>%
    group_by(outcome, subsample, ev_yr) 
# Add new columns for color and shape mapping
chainindepcontrols$color <- factor(chainindepcontrols$subsample)
chainindepcontrols$shape <- factor(chainindepcontrols$subsample)
chainindepcontrolspre <- filter(chainindepcontrols, ev_yr < 0) %>%
    group_by(outcome, subsample, ev_yr) 
chainindepcontrolspost <- filter(chainindepcontrols, ev_yr >= 0) %>%
    group_by(outcome, subsample, ev_yr) 

# Create the chain/independent plot
chainindepcontrols_plot <- ggplot(data = chainindepcontrols, aes(x = ev_yr, y = b, color = color, shape = shape)) +
    geom_linerange(aes(ymin = ll, ymax = ul),
                   alpha = 0.7,
                   size = 0.8) +
    geom_point(size = 3) +
    geom_vline(xintercept = 0, lty = 3, colour = "red") + 
    geom_hline(yintercept = 0, lty = 2) +
    scale_color_manual(values = colors, name = "Chain status") +
    scale_shape_manual(values = shapes, name = "Chain status") +
    geom_line(data = chainindepcontrolspre, aes(y = avgb, x = ev_yr, color = color)) + 
    geom_ribbon(data = chainindepcontrolspre, aes(ymin = avgll, ymax = avgul, color = color),
                   alpha = 0.1) +
    geom_line(data = chainindepcontrolspost, aes(y = avgb, x = ev_yr, color = color)) + 
    geom_ribbon(data = chainindepcontrolspost, aes(ymin = avgll, ymax = avgul, color = color),
                   alpha = 0.1) +
    labs(y = "ATT estimate", x = "Event Year", 
    caption = paste0("ATT estimate for the chain subsample in the post period is ",
                as.character(round(filter(chainindepcontrolspost, subsample == "Chain")$avgb[1] , digits = 3)),
                " with a 95% confidence\ninterval of (",
                as.character(round(filter(chainindepcontrolspost, subsample == "Chain")$avgll[1], digits = 3)),
                ", ",
                as.character(round(filter(chainindepcontrolspost, subsample == "Chain")$avgul[1], digits = 3)),
                ").\nATT estimate for the independent subsample in the post period is ",
                as.character(round(filter(chainindepcontrolspost, subsample == "Independent")$avgb[1], digits = 3)),
                " with a 95% confidence\ninterval of (",
                as.character(round(filter(chainindepcontrolspost, subsample == "Independent")$avgll[1], digits = 3)),
                ", ",
                as.character(round(filter(chainindepcontrolspost, subsample == "Independent")$avgul[1], digits = 3)),
                ").\nCounty level controls - median income, population, population < 5, share Hispanic,\nshare Black, share earning < FPL, share earning <180% FPL, and share receiving\nSNAP or cash benefits.")) +
    theme_light() +
    scale_x_continuous(breaks = seq(-4, 4, 2)) +
    theme(plot.caption = element_text(hjust = 0), legend.position="none") + 
    facet_wrap(~ subsample, ncol = 2)
  
ggsave("csdid_tip_auth_chainindepcontrols.png", plot = chainindepcontrols_plot, width = 6.5, height = 4, units = "in")

# make processor graphs
# load in implementation years to use when figuring out implementation order
ewic <- read_dta("./data/cleaned/ewic_rollout.dta") %>%
  select(c(state, ev_year, ct_fips)) %>%
  distinct()

tip <- read_dta("./data/cleaned/tip_auth_sq.dta") %>%
    select(c(state, ct_fips)) %>%
    distinct() %>%
  rename(st = state)

process <- left_join(tip, ewic, by = "ct_fips") %>%
  filter(!((st == "KY" & state == "Tennessee") | (st == "NH" & state == "Maine"))) %>%
    mutate(processor = case_when(st %in% c("VA", "CT", "NH", "NY", "MI", "OK", "RI", "AL", "GA", "MS", "SC", "TN", "IN") ~ "Conduent",
                            st %in% c("NM", "TX", "ME", "MD", "NJ", "PA", "NC", "IL", "AR", "LA", "UT", "MT") ~ "Solutran", #note that solutran implemented wyoming in 2002 before new mexico or texas
                            st %in% c("MA", "WV", "FL", "KY", "WI", "OR", "DE", "DC", "IA", "MI", "MN", "AZ", "CO", "KS", "MO", "NE") ~ "FIS CDP")) %>%
    group_by(st) %>%
    mutate(ev_year = min(ev_year)) %>%
    mutate(avg_ev_year = mean(ev_year)) %>%
    select(c(st, ev_year, processor, avg_ev_year)) %>%
    rename(state = st) %>%
    distinct()

simple <- read_dta("./analysis/output/csdid_tip_simple.dta") 

processor_all <- filter(simple, state != "all" & subsample == "all" & outcome == "auth" & controls == "no") %>%
    left_join(process, by = "state") %>%
    group_by(processor) %>%
    arrange(processor, ev_year, avg_ev_year) %>%
    mutate(n = row_number()) %>%
    filter(!is.na(processor) & !is.na(se))

# Define the colors for each chain status
colors <- brewer.pal(3, "Set2")

# Define the shapes for each chain status
shapes <- c(16, 17, 18)

# Add new columns for color and shape mapping
processor_all$color <- factor(processor_all$processor)
processor_all$shape <- factor(processor_all$processor)

# Create the chain/independent plot
processor_all_plot <- ggplot(data = processor_all, aes(x = n, y = b, color = color, shape = shape)) +
    geom_linerange(aes(ymin = ll, ymax = ul),
                   alpha = 0.7,
                   size = 0.8) +
    geom_point(size = 3) +
    geom_smooth(method = lm, se = FALSE) +
    stat_poly_eq(formula = y ~ x,
                output.type = "numeric",
                mapping = aes(label = sprintf("Slope: %s",
                            c(formatC(round(stat(coef.ls)[[1]][[2, "Estimate"]], digits = 3)),
                              formatC(round(stat(coef.ls)[[2]][[2, "Estimate"]], digits = 3)),
                              formatC(round(stat(coef.ls)[[3]][[2, "Estimate"]], digits = 3))))),
                label.x = 3.5,
                label.y = -0.05) +
    geom_hline(yintercept = 0, lty = 2) +
    scale_color_manual(values = colors, name = "Processor") +
    scale_shape_manual(values = shapes, name = "Processor") +
    labs(y = "ATT estimate", x = "Processor order", 
    caption = "Graph shows ATT estimates for individual states, grouped by the processor that implemented\nWIC EBT in the state.\nThe ordering on the x-axis represents the order in which processors worked with the states,\ndetermined by the year of first EBT implementation in the state.\nWe assume that earlier implementation results from being earlier in working with the processor.") +
    theme_light() +
    theme(plot.caption = element_text(hjust = 0), legend.position="none") + 
    facet_wrap(~ processor, ncol = 2)

ggsave("csdid_tip_auth_processor.png", plot = processor_all_plot, width = 6.5, height = 4, units = "in")

processorcontrols_all <- filter(simple, state != "all" & subsample == "all" & outcome == "auth" & controls == "yes") %>%
    left_join(process, by = "state") %>%
    group_by(processor) %>%
    arrange(processor, ev_year, avg_ev_year) %>%
    mutate(n = row_number()) %>%
    filter(!is.na(processor) & !is.na(se))

# Define the colors for each chain status
colors <- brewer.pal(3, "Set2")

# Define the shapes for each chain status
shapes <- c(16, 17, 18)

# Add new columns for color and shape mapping
processorcontrols_all$color <- factor(processorcontrols_all$processor)
processorcontrols_all$shape <- factor(processorcontrols_all$processor)
# Create the chain/independent plot
processorcontrols_all_plot <- ggplot(data = processorcontrols_all, aes(x = n, y = b, color = color, shape = shape)) +
    geom_linerange(aes(ymin = ll, ymax = ul),
                   alpha = 0.7,
                   size = 0.8) +
    geom_point(size = 3) +
    geom_smooth(method = lm, se = FALSE) +
    stat_poly_eq(formula = y ~ x,
                output.type = "numeric",
                mapping = aes(label = sprintf("Slope: %s",
                            c(formatC(round(stat(coef.ls)[[1]][[2, "Estimate"]], digits = 3)),
                              formatC(round(stat(coef.ls)[[2]][[2, "Estimate"]], digits = 3)),
                              formatC(round(stat(coef.ls)[[3]][[2, "Estimate"]], digits = 3))))),
                label.x = 3.5,
                label.y = -0.05) +
    geom_hline(yintercept = 0, lty = 2) +
    scale_color_manual(values = colors, name = "Processor") +
    scale_shape_manual(values = shapes, name = "Processor") +
    labs(y = "ATT estimate", x = "Processor order", 
    caption = "Graph shows ATT estimates for individual states, grouped by the processor that implemented\nWIC EBT in the state.\nThe ordering on the x-axis represents the order in which processors worked with the states,\ndetermined by the year of first EBT implementation in the state.\nWe assume that earlier implementation results from being earlier in working with the processor.\nCounty level controls - median income, population, population < 5, share Hispanic,\nshare Black, share earning < FPL, share earning <180% FPL, and share receiving\nSNAP or cash benefits.") +
    theme_light() +
    theme(plot.caption = element_text(hjust = 0), legend.position="none") + 
    facet_wrap(~ processor, ncol = 2)

ggsave("csdid_tip_auth_processor_controls.png", plot = processorcontrols_all_plot, width = 6.5, height = 4, units = "in")

processor_chain <-  filter(simple, state != "all" & subsample == "chain" & outcome == "auth" & controls == "no") %>%
    left_join(process, by = "state") %>%
    group_by(processor) %>%
    arrange(processor, ev_year, avg_ev_year) %>%
    mutate(n = row_number()) %>%
  filter(!is.na(processor) & !is.na(se))

# Add new columns for color and shape mapping
processor_chain$color <- factor(processor_chain$processor)
processor_chain$shape <- factor(processor_chain$processor)
# Create the chain/independent plot
processor_chain_plot <- ggplot(data = processor_chain, aes(x = n, y = b, color = color, shape = shape)) +
    geom_linerange(aes(ymin = ll, ymax = ul),
                   alpha = 0.7,
                   size = 0.8) +
    geom_point(size = 3) +
    geom_smooth(method = lm, se = FALSE) +
    stat_poly_eq(formula = y ~ x,
                output.type = "numeric",
                mapping = aes(label = sprintf("Slope: %s",
                            c(formatC(round(stat(coef.ls)[[1]][[2, "Estimate"]], digits = 3)),
                              formatC(round(stat(coef.ls)[[2]][[2, "Estimate"]], digits = 3)),
                              formatC(round(stat(coef.ls)[[3]][[2, "Estimate"]], digits = 3))))),
                label.x = 3.5,
                label.y = -0.05) +
    geom_hline(yintercept = 0, lty = 2) +
    scale_color_manual(values = colors, name = "Processor") +
    scale_shape_manual(values = shapes, name = "Processor") +
    labs(y = "ATT estimate, chain stores", x = "Processor order", 
    caption = "Graph shows ATT estimates for individual states' chain stores, grouped by the processor that\nimplemented WIC EBT in the state.\nThe ordering on the x-axis represents the order in which processors worked with the states,\ndetermined by the year of first EBT implementation in the state.\nWe assume that earlier implementation results from being earlier in working with the processor.") +
    theme_light() +
    theme(plot.caption = element_text(hjust = 0), legend.position="none") + 
    facet_wrap(~ processor, ncol = 1)

ggsave("csdid_tip_auth_processor_chain.png", plot = processor_chain_plot, width = 6.5, units = "in")

processor_chain_controls <-  filter(simple, state != "all" & subsample == "chain" & outcome == "auth" & controls == "yes") %>%
    left_join(process, by = "state") %>%
    group_by(processor) %>%
    arrange(processor, ev_year, avg_ev_year) %>%
    mutate(n = row_number()) %>%
  filter(!is.na(processor) & !is.na(se))

# Add new columns for color and shape mapping
processor_chain_controls$color <- factor(processor_chain_controls$processor)
processor_chain_controls$shape <- factor(processor_chain_controls$processor)
# Create the chain/independent plot
processor_chain_controls_plot <- ggplot(data = processor_chain_controls, aes(x = n, y = b, color = color, shape = shape)) +
    geom_linerange(aes(ymin = ll, ymax = ul),
                   alpha = 0.7,
                   size = 0.8) +
    geom_point(size = 3) +
    geom_smooth(method = lm, se = FALSE) +
    stat_poly_eq(formula = y ~ x,
                output.type = "numeric",
                mapping = aes(label = sprintf("Slope: %s",
                            c(formatC(round(stat(coef.ls)[[1]][[2, "Estimate"]], digits = 3)),
                              formatC(round(stat(coef.ls)[[2]][[2, "Estimate"]], digits = 3)),
                              formatC(round(stat(coef.ls)[[3]][[2, "Estimate"]], digits = 3))))),
                label.x = 3.5,
                label.y = -0.05) +
    geom_hline(yintercept = 0, lty = 2) +
    scale_color_manual(values = colors, name = "Processor") +
    scale_shape_manual(values = shapes, name = "Processor") +
    labs(y = "ATT estimate, chain stores", x = "Processor order", 
    caption = "Graph shows ATT estimates for individual states' chain stores, grouped by the processor that\nimplemented WIC EBT in the state.\nThe ordering on the x-axis represents the order in which processors worked with the states,\ndetermined by the year of first EBT implementation in the state.\nWe assume that earlier implementation results from being earlier in working with the processor.\nCounty level controls - median income, population, population < 5, share Hispanic,\nshare Black, share earning < FPL, share earning <180% FPL, and share receiving\nSNAP or cash benefits.") +
    theme_light() +
    theme(plot.caption = element_text(hjust = 0), legend.position="none") + 
    facet_wrap(~ processor, ncol = 1)

ggsave("csdid_tip_auth_processor_chain_controls.png", plot = processor_chain_controls_plot, width = 6.5, units = "in")

processor_indep <-  filter(simple, state != "all" & subsample == "indep" & outcome == "auth" & controls == "no") %>%
    left_join(process, by = "state") %>%
    group_by(processor) %>%
    arrange(processor, ev_year, avg_ev_year) %>%
    mutate(n = row_number()) %>%
  filter(!is.na(processor) & !is.na(se))

# Add new columns for color and shape mapping
processor_indep$color <- factor(processor_indep$processor)
processor_indep$shape <- factor(processor_indep$processor)
# Create the chain/independent plot
processor_indep_plot <- ggplot(data = processor_indep, aes(x = n, y = b, color = color, shape = shape)) +
    geom_linerange(aes(ymin = ll, ymax = ul),
                   alpha = 0.7,
                   size = 0.8) +
    geom_point(size = 3) +
    geom_hline(yintercept = 0, lty = 2) +
    geom_smooth(method = lm, se = FALSE) +
    stat_poly_eq(formula = y ~ x,
                output.type = "numeric",
                mapping = aes(label = sprintf("Slope: %s",
                            c(formatC(round(stat(coef.ls)[[1]][[2, "Estimate"]], digits = 3)),
                              formatC(round(stat(coef.ls)[[2]][[2, "Estimate"]], digits = 3)),
                              formatC(round(stat(coef.ls)[[3]][[2, "Estimate"]], digits = 3))))),
                label.x = 3.5,
                label.y = -0.05) +
    scale_color_manual(values = colors, name = "Processor") +
    scale_shape_manual(values = shapes, name = "Processor") +
    labs(y = "ATT estimate, independent stores", x = "Processor order", 
    caption = "Graph shows ATT estimates for individual states' independent stores, grouped by the\nsprocessor that implemented WIC EBT in the state.\nThe ordering on the x-axis represents the order in which processors worked with the states,\ndetermined by the year of first EBT implementation in the state.\nWe assume that earlier implementation results from being earlier in working with the processor.") +
    theme_light() +
    theme(plot.caption = element_text(hjust = 0), legend.position="none") +  
    facet_wrap(~ processor, ncol = 1)

ggsave("csdid_tip_auth_processor_indep.png", plot = processor_indep_plot, width = 6.5, units = "in")

processor_indep_controls <-  filter(simple, state != "all" & subsample == "indep" & outcome == "auth" & controls == "yes") %>%
    left_join(process, by = "state") %>%
    group_by(processor) %>%
    arrange(processor, ev_year, avg_ev_year) %>%
    mutate(n = row_number()) %>%
  filter(!is.na(processor) & !is.na(se))

# Add new columns for color and shape mapping
processor_indep_controls$color <- factor(processor_indep_controls$processor)
processor_indep_controls$shape <- factor(processor_indep_controls$processor)
# Create the chain/independent plot
processor_indep_controls_plot <- ggplot(data = processor_indep_controls, aes(x = n, y = b, color = color, shape = shape)) +
    geom_linerange(aes(ymin = ll, ymax = ul),
                   alpha = 0.7,
                   size = 0.8) +
    geom_point(size = 3) +
    geom_smooth(method = lm, se = FALSE) +
    stat_poly_eq(formula = y ~ x,
                output.type = "numeric",
                mapping = aes(label = sprintf("Slope: %s",
                            c(formatC(round(stat(coef.ls)[[1]][[2, "Estimate"]], digits = 3)),
                              formatC(round(stat(coef.ls)[[2]][[2, "Estimate"]], digits = 3)),
                              formatC(round(stat(coef.ls)[[3]][[2, "Estimate"]], digits = 3))))),
                label.x = 3.5,
                label.y = -0.05) +
    geom_hline(yintercept = 0, lty = 2) +
    scale_color_manual(values = colors, name = "Processor") +
    scale_shape_manual(values = shapes, name = "Processor") +
    labs(y = "ATT estimate, independent stores", x = "Processor order", 
    caption = "Graph shows ATT estimates for individual states' independent stores, grouped by the\nsprocessor that implemented WIC EBT in the state.\nThe ordering on the x-axis represents the order in which processors worked with the states,\ndetermined by the year of first EBT implementation in the state.\nWe assume that earlier implementation results from being earlier in working with the processor.\nCounty level controls - median income, population, population < 5, share Hispanic,\nshare Black, share earning < FPL, share earning <180% FPL, and share receiving\nSNAP or cash benefits.") +
    theme_light() +
    theme(plot.caption = element_text(hjust = 0), legend.position="none") +  
    facet_wrap(~ processor, ncol = 1)

ggsave("csdid_tip_auth_processor_indep_controls.png", plot = processor_indep_controls_plot, width = 6.5, units = "in")

### make treatment group specific graphs
group <- read_dta("./analysis/output/csdid_tip_group.dta") 

picture <- filter(group, state == "all") %>%
  group_by(outcome, subsample, gr_yr, controls) %>%
  mutate(subsample = gsub("all","All", subsample)) %>%
  mutate(subsample = gsub("chain","Chain", subsample)) %>%
  mutate(subsample = gsub("indep","Independent", subsample)) 

avg <- filter(picture, gr_yr == "Average") %>%
    ungroup() %>%
    select(!c(gr_yr)) %>%
    mutate(avgb = b, avgll = ll, avgul = ul) %>%
    select(c(outcome, subsample, state, controls, avgb, avgll, avgul))

picture <- filter(picture, gr_yr != "Average") %>%
  mutate(gr_yr = as.numeric(gr_yr))
picture <- left_join(picture, avg, by = c("outcome", "subsample", "state", "controls"))

rm(avg)

# make three pictures:
#one that just shows all
#one that compares chain and indep
#one that compares chain and indep with different timings

setwd("./analysis/output/graphs")

all <- filter(picture, subsample == "All" & controls == "no")

group_plot <- ggplot(all, aes(x = gr_yr, y = b)) +
    geom_point(size = 3) +
    geom_ribbon(aes(ymin = all$avgll[1], ymax = all$avgul[1]),
              alpha = 0.25,
              color = "gray90") +
    geom_linerange(aes(ymin = ll, ymax = ul),
                alpha = 0.7,
                size = 0.8) +
    geom_hline(yintercept = 0, lty = 2) +
    geom_hline(yintercept = all$avgb[1], color = "gray70") + 
    labs(y = "ATT estimate",
    x = "FY WIC EBT First Implemented",
    caption = paste0("ATT estimate across all groups is ",
                     as.character(round(all$avgb[1], digits = 3)),
                     " with a 95% confidence interval of (",
                     as.character(round(all$avgll[1], digits = 3)),
                     ", ",
                     as.character(round(all$avgul[1], digits = 3)),
                     ").")) +
    theme_light() + 
    scale_x_continuous(breaks = seq(2006, 2016, 2)) +
    theme(plot.caption = element_text(hjust = 0))
  
ggsave("csdid_tip_auth_all_group.png", plot = group_plot, width = 6.5, height = 4, units = "in")

group_plot <- ggplot(all, aes(x = gr_yr, y = b)) +
    geom_point(size = 3) +
    geom_ribbon(aes(ymin = all$avgll[1], ymax = all$avgul[1]),
              alpha = 0.25,
              color = "gray90") +
    geom_linerange(aes(ymin = ll, ymax = ul),
                alpha = 0.7,
                size = 0.8) +
    geom_smooth(method = lm, se = FALSE) +
    stat_poly_eq(formula = y ~ x,
                output.type = "numeric",
                mapping = aes(label = sprintf("Slope: %s",
                            c(formatC(round(stat(coef.ls)[[1]][[2, "Estimate"]], digits = 3))))),
                label.x = 3.5,
                label.y = -0.05) +
    geom_hline(yintercept = 0, lty = 2) +
    geom_hline(yintercept = all$avgb[1], color = "gray70") + 
    labs(y = "ATT estimate",
    x = "FY WIC EBT First Implemented",
    caption = paste0("ATT estimate across all groups is ",
                     as.character(round(all$avgb[1], digits = 3)),
                     " with a 95% confidence interval of (",
                     as.character(round(all$avgll[1], digits = 3)),
                     ", ",
                     as.character(round(all$avgul[1], digits = 3)),
                     ").")) +
    theme_light() + 
    scale_x_continuous(breaks = seq(2006, 2016, 2)) +
    theme(plot.caption = element_text(hjust = 0))
  
ggsave("csdid_tip_auth_all_group_fitline.png", plot = group_plot, width = 6.5, height = 4, units = "in")

# Define the colors for each chain status
colors <- brewer.pal(3, "Set2")[1:2]

# Define the shapes for each chain status
shapes <- c(16, 17)

chainindep <- filter(picture, controls == "no" & (subsample == "Chain" | subsample == "Independent")) %>%
    group_by(subsample) 
# Add new columns for color and shape mapping
chainindep$color <- factor(chainindep$subsample)
chainindep$shape <- factor(chainindep$subsample)

# Create the chain/independent plot
chainindep_plot <- ggplot(data = chainindep, aes(x = gr_yr, y = b, color = color, shape = shape)) +
    geom_linerange(aes(ymin = ll, ymax = ul),
                   alpha = 0.7,
                   size = 0.8) +
    geom_point(size = 3) +
    geom_hline(yintercept = 0, lty = 2) +
    scale_color_manual(values = colors, name = "Chain status") +
    scale_shape_manual(values = shapes, name = "Chain status") +
    geom_line(data = chainindep, aes(y = avgb, x = gr_yr, color = color)) + 
    geom_ribbon(data = chainindep, aes(ymin = avgll, ymax = avgul, color = color),
                   alpha = 0.1) +
    labs(y = "ATT estimate", x = "FY WIC EBT First Implemented", 
    caption = paste0("ATT estimate for the chain subsample is ",
                as.character(round(filter(chainindep, subsample == "Chain")$avgb[1] , digits = 3)),
                " with a 95% confidence interval of (",
                as.character(round(filter(chainindep, subsample == "Chain")$avgll[1], digits = 3)),
                ", ",
                as.character(round(filter(chainindep, subsample == "Chain")$avgul[1], digits = 3)),
                ").\nATT estimate for the independent subsample is ",
                as.character(round(filter(chainindep, subsample == "Independent")$avgb[1], digits = 3)),
                " with a 95% confidence interval of\n(",
                as.character(round(filter(chainindep, subsample == "Independent")$avgll[1], digits = 3)),
                ", ",
                as.character(round(filter(chainindep, subsample == "Independent")$avgul[1], digits = 3)),
                ").")) +
    theme_light() +
    theme(plot.caption = element_text(hjust = 0), legend.position="none") + 
    scale_x_continuous(breaks = seq(2006, 2016, 2)) +
    facet_wrap(~ subsample, ncol = 2)
  
ggsave("csdid_tip_auth_chainindep_group.png", plot = chainindep_plot, width = 6.5, height = 4, units = "in")


# Create the chain/independent plot with fitline
chainindep_plot <- ggplot(data = chainindep, aes(x = gr_yr, y = b, color = color, shape = shape)) +
    geom_linerange(aes(ymin = ll, ymax = ul),
                   alpha = 0.7,
                   size = 0.8) +
    geom_point(size = 3) +
    geom_smooth(method = lm, se = FALSE) +
    stat_poly_eq(formula = y ~ x,
                output.type = "numeric",
                mapping = aes(label = sprintf("Slope: %s",
                            c(formatC(round(stat(coef.ls)[[1]][[2, "Estimate"]], digits = 3)),
                              formatC(round(stat(coef.ls)[[2]][[2, "Estimate"]], digits = 3))))),
                label.x = 3.5,
                label.y = -0.05) +
    geom_hline(yintercept = 0, lty = 2) +
    scale_color_manual(values = colors, name = "Chain status") +
    scale_shape_manual(values = shapes, name = "Chain status") +
    geom_line(data = chainindep, aes(y = avgb, x = gr_yr, color = color)) + 
    geom_ribbon(data = chainindep, aes(ymin = avgll, ymax = avgul, color = color),
                   alpha = 0.1) +
    labs(y = "ATT estimate", x = "FY WIC EBT First Implemented", 
    caption = paste0("ATT estimate for the chain subsample is ",
                as.character(round(filter(chainindep, subsample == "Chain")$avgb[1] , digits = 3)),
                " with a 95% confidence interval of (",
                as.character(round(filter(chainindep, subsample == "Chain")$avgll[1], digits = 3)),
                ", ",
                as.character(round(filter(chainindep, subsample == "Chain")$avgul[1], digits = 3)),
                ").\nATT estimate for the independent subsample is ",
                as.character(round(filter(chainindep, subsample == "Independent")$avgb[1], digits = 3)),
                " with a 95% confidence interval of\n(",
                as.character(round(filter(chainindep, subsample == "Independent")$avgll[1], digits = 3)),
                ", ",
                as.character(round(filter(chainindep, subsample == "Independent")$avgul[1], digits = 3)),
                ").")) +
    theme_light() +
    theme(plot.caption = element_text(hjust = 0), legend.position="none") + 
    scale_x_continuous(breaks = seq(2006, 2016, 2)) +
    facet_wrap(~ subsample, ncol = 2)
  
ggsave("csdid_tip_auth_chainindep_group_fitline.png", plot = chainindep_plot, width = 6.5, height = 4, units = "in")


all <- filter(picture, subsample == "All" & controls == "yes")

group_plot <- ggplot(all, aes(x = gr_yr, y = b)) +
    geom_point(size = 3) +
    geom_ribbon(aes(ymin = all$avgll[1], ymax = all$avgul[1]),
              alpha = 0.25,
              color = "gray90") +
    geom_linerange(aes(ymin = ll, ymax = ul),
                alpha = 0.7,
                size = 0.8) +
    geom_hline(yintercept = 0, lty = 2) +
    geom_hline(yintercept = all$avgb[1], color = "gray70") + 
    labs(y = "ATT estimate",
    x = "FY WIC EBT First Implemented",
    caption = paste0("ATT estimate across all groups is ",
                     as.character(round(all$avgb[1], digits = 3)),
                     " with a 95% confidence interval of (",
                     as.character(round(all$avgll[1], digits = 3)),
                     ", ",
                     as.character(round(all$avgul[1], digits = 3)),
                     ").\nCounty level controls - median income, population, population < 5, share Hispanic,\nshare Black, share earning < FPL, share earning <180% FPL, and share receiving\nSNAP or cash benefits.")) +
    theme_light() + 
    scale_x_continuous(breaks = seq(2006, 2016, 2)) +
    theme(plot.caption = element_text(hjust = 0))
  
ggsave("csdid_tip_auth_all_group_controls.png", plot = group_plot, width = 6.5, height = 4, units = "in")

# Define the colors for each chain status
colors <- brewer.pal(3, "Set2")[1:2]

# Define the shapes for each chain status
shapes <- c(16, 17)

chainindep <- filter(picture, controls == "yes" & (subsample == "Chain" | subsample == "Independent")) %>%
    group_by(subsample) 
# Add new columns for color and shape mapping
chainindep$color <- factor(chainindep$subsample)
chainindep$shape <- factor(chainindep$subsample)

# Create the chain/independent plot
chainindep_plot <- ggplot(data = chainindep, aes(x = gr_yr, y = b, color = color, shape = shape)) +
    geom_linerange(aes(ymin = ll, ymax = ul),
                   alpha = 0.7,
                   size = 0.8) +
    geom_point(size = 3) +
    geom_hline(yintercept = 0, lty = 2) +
    scale_color_manual(values = colors, name = "Chain status") +
    scale_shape_manual(values = shapes, name = "Chain status") +
    geom_line(data = chainindep, aes(y = avgb, x = gr_yr, color = color)) + 
    geom_ribbon(data = chainindep, aes(ymin = avgll, ymax = avgul, color = color),
                   alpha = 0.1) +
    labs(y = "ATT estimate", x = "FY WIC EBT First Implemented", 
    caption = paste0("ATT estimate for the chain subsample is ",
                as.character(round(filter(chainindep, subsample == "Chain")$avgb[1] , digits = 3)),
                " with a 95% confidence interval of (",
                as.character(round(filter(chainindep, subsample == "Chain")$avgll[1], digits = 3)),
                ", ",
                as.character(round(filter(chainindep, subsample == "Chain")$avgul[1], digits = 3)),
                ").\nATT estimate for the independent subsample is ",
                as.character(round(filter(chainindep, subsample == "Independent")$avgb[1], digits = 3)),
                " with a 95% confidence interval of\n(",
                as.character(round(filter(chainindep, subsample == "Independent")$avgll[1], digits = 3)),
                ", ",
                as.character(round(filter(chainindep, subsample == "Independent")$avgul[1], digits = 3)),
                ").\nCounty level controls - median income, population, population < 5, share Hispanic,\nshare Black, share earning < FPL, share earning <180% FPL, and share receiving\nSNAP or cash benefits.")) +
    theme_light() +
    theme(plot.caption = element_text(hjust = 0), legend.position="none") + 
    scale_x_continuous(breaks = seq(2006, 2016, 2)) +
    facet_wrap(~ subsample, ncol = 2)
  
ggsave("csdid_tip_auth_chainindep_group_controls.png", plot = chainindep_plot, width = 6.5, height = 4, units = "in")

#### switch to zip level analyses (of WIC redemptions )
event <- read_dta("./analysis/output/csdid_zip_event.dta") 

picture <- event %>%
  group_by(outcome, subsample, ev_yr) %>%
  mutate(subsample = gsub("all","All", subsample)) %>%
  mutate(subsample = gsub("highpov","High poverty ZIPs", subsample)) 

# make one picture:
#one that just shows all

setwd("./analysis/output/graphs")

all <- filter(picture, subsample == "All")
allpre <- filter(picture, subsample == "All", ev_yr < 0)
allpost <- filter(picture, subsample == "All", ev_yr >= 0)

all_plot <- ggplot(all, aes(x = ev_yr, y = b)) +
    geom_linerange(aes(ymin = ll, ymax = ul),
                   alpha = 0.7,
                   size = 0.8) +
    geom_point(size = 3) +
    geom_vline(xintercept = 0, lty = 3, colour = "red") + 
    geom_hline(yintercept = 0, lty = 2) +
    geom_line(data = allpre, aes(x = ev_yr, y = avgb), colour = "gray70") + 
    geom_ribbon(data = allpre, aes(ymin = avgll, ymax = avgul),
                   alpha = 0.25,
                   colour = "gray90") +
    geom_line(data = allpost, aes(x = ev_yr, y = avgb), colour = "gray70") + 
    geom_ribbon(data = allpost, aes(ymin = avgll, ymax = avgul),
                   alpha = 0.25,
                   colour = "gray90") +
    labs(y = "ATT estimate",
    x = "Event Year",
    caption = paste0("ATT estimate for the post period shown here is ",
                     as.character(round(allpost$avgb[1], digits = 3)),
                     " with a 95% confidence interval of (",
                     as.character(round(allpost$avgll[1], digits = 3)),
                     ", ",
                     as.character(round(allpost$avgul[1], digits = 3)),
                     ").")) +
    theme_light() + 
    scale_x_continuous(breaks = seq(-4, 4, 2)) +
  theme(plot.caption = element_text(hjust = 0), legend.position="none") 
  
ggsave("csdid_zip_wicredhk_all.png", plot = all_plot)

# pictures for all for TIP for ZIPs only
setwd("~/Shared/pc-ewic/analysis/output/graphs/")
event <- read_dta("./analysis/output/csdid_tip_4zips_event.dta") 

picture <- filter(event, state == "all") %>%
  group_by(outcome, subsample, controls, ev_yr) %>%
  mutate(subsample = gsub("all","All", subsample)) %>%
  mutate(subsample = gsub("chain","Chain", subsample)) %>%
  mutate(subsample = gsub("indep","Independent", subsample)) %>%
  mutate(controls = gsub("no", "No", controls)) %>%
  mutate(controls = gsub("yes", "Yes", controls))

all <- filter(picture, subsample == "All" & controls == "No")
allpre <- filter(picture, subsample == "All" & controls == "No", ev_yr < 0)
allpost <- filter(picture, subsample == "All" & controls == "No", ev_yr >= 0)

all_plot <- ggplot(all, aes(x = ev_yr, y = b)) +
    geom_linerange(aes(ymin = ll, ymax = ul),
                   alpha = 0.7,
                   size = 0.8) +
    geom_point(size = 3) +
    geom_vline(xintercept = 0, lty = 3, colour = "red") + 
    geom_hline(yintercept = 0, lty = 2) +
    geom_line(data = allpre, aes(x = ev_yr, y = avgb), colour = "gray70") + 
    geom_ribbon(data = allpre, aes(ymin = avgll, ymax = avgul),
                   alpha = 0.25,
                   colour = "gray90") +
    geom_line(data = allpost, aes(x = ev_yr, y = avgb), colour = "gray70") + 
    geom_ribbon(data = allpost, aes(ymin = avgll, ymax = avgul),
                   alpha = 0.25,
                   colour = "gray90") +
    labs(y = "ATT estimate",
    x = "Event Year",
    caption = paste0("ATT estimate for the post period shown here is ",
                     as.character(round(allpost$avgb[1], digits = 3)),
                     " with a 95% confidence interval of\n(",
                     as.character(round(allpost$avgll[1], digits = 3)),
                     ", ",
                     as.character(round(allpost$avgul[1], digits = 3)),
                     ").")) +
    theme_light() +
    scale_x_continuous(breaks = seq(-4, 4, 2)) +
    theme(plot.caption = element_text(hjust = 0))
  
ggsave("csdid_tip_4zips_auth_all.png", plot = all_plot, width = 6.5, height = 4, units = "in")

# Define the colors for each chain status
colors <- brewer.pal(3, "Set2")[1:2]

# Define the shapes for each chain status
shapes <- c(16, 17)

chainindep <- filter(picture, (subsample == "Chain" | subsample == "Independent") & controls == "No") %>%
    group_by(outcome, subsample, ev_yr) 
# Add new columns for color and shape mapping
chainindep$color <- factor(chainindep$subsample)
chainindep$shape <- factor(chainindep$subsample)
chainindeppre <- filter(chainindep, ev_yr < 0) %>%
    group_by(outcome, subsample, ev_yr) 
chainindeppost <- filter(chainindep, ev_yr >= 0) %>%
    group_by(outcome, subsample, ev_yr) 

# Create the chain/independent plot
chainindep_plot <- ggplot(data = chainindep, aes(x = ev_yr, y = b, color = color, shape = shape)) +
    geom_linerange(aes(ymin = ll, ymax = ul),
                   alpha = 0.7,
                   size = 0.8) +
    geom_point(size = 3) +
    geom_vline(xintercept = 0, lty = 3, colour = "red") + 
    geom_hline(yintercept = 0, lty = 2) +
    scale_color_manual(values = colors, name = "Chain status") +
    scale_shape_manual(values = shapes, name = "Chain status") +
    geom_line(data = chainindeppre, aes(y = avgb, x = ev_yr, color = color)) + 
    geom_ribbon(data = chainindeppre, aes(ymin = avgll, ymax = avgul, color = color),
                   alpha = 0.1) +
    geom_line(data = chainindeppost, aes(y = avgb, x = ev_yr, color = color)) + 
    geom_ribbon(data = chainindeppost, aes(ymin = avgll, ymax = avgul, color = color),
                   alpha = 0.1) +
    labs(y = "ATT estimate", x = "Event Year", 
    caption = paste0("ATT estimate for the chain subsample in the post period is ",
                as.character(round(filter(chainindeppost, subsample == "Chain")$avgb[1] , digits = 3)),
                " with a 95% confidence\ninterval of (",
                as.character(round(filter(chainindeppost, subsample == "Chain")$avgll[1], digits = 3)),
                ", ",
                as.character(round(filter(chainindeppost, subsample == "Chain")$avgul[1], digits = 3)),
                ").\nATT estimate for the independent subsample in the post period is ",
                as.character(round(filter(chainindeppost, subsample == "Independent")$avgb[1], digits = 3)),
                " with a 95% confidence\ninterval of (",
                as.character(round(filter(chainindeppost, subsample == "Independent")$avgll[1], digits = 3)),
                ", ",
                as.character(round(filter(chainindeppost, subsample == "Independent")$avgul[1], digits = 3)),
                ").")) +
    theme_light() +
    scale_x_continuous(breaks = seq(-4, 4, 2)) +
    theme(plot.caption = element_text(hjust = 0), legend.position="none") + 
    facet_wrap(~ subsample, ncol = 2)
  
ggsave("csdid_tip_4zips_auth_chainindep.png", plot = chainindep_plot, width = 6.5, height = 4, units = "in")

### make treatment group specific graphs
group <- read_dta("./analysis/output/csdid_tip_4zips_group.dta") 

picture <- filter(group, state == "all") %>%
  group_by(outcome, subsample, gr_yr, controls) %>%
  mutate(subsample = gsub("all","All", subsample)) %>%
  mutate(subsample = gsub("chain","Chain", subsample)) %>%
  mutate(subsample = gsub("indep","Independent", subsample)) 

avg <- filter(picture, gr_yr == "Average") %>%
    ungroup() %>%
    select(!c(gr_yr)) %>%
    mutate(avgb = b, avgll = ll, avgul = ul) %>%
    select(c(outcome, subsample, state, controls, avgb, avgll, avgul))

picture <- filter(picture, gr_yr != "Average") %>%
  mutate(gr_yr = as.numeric(gr_yr))
picture <- left_join(picture, avg, by = c("outcome", "subsample", "state", "controls"))

rm(avg)

# make three pictures:
#one that just shows all
#one that compares chain and indep
#one that compares chain and indep with different timings

setwd("./analysis/output/graphs")

all <- filter(picture, subsample == "All" & controls == "no")

group_plot <- ggplot(all, aes(x = gr_yr, y = b)) +
    geom_point(size = 3) +
    geom_ribbon(aes(ymin = all$avgll[1], ymax = all$avgul[1]),
              alpha = 0.25,
              color = "gray90") +
    geom_linerange(aes(ymin = ll, ymax = ul),
                alpha = 0.7,
                size = 0.8) +
    geom_hline(yintercept = 0, lty = 2) +
    geom_hline(yintercept = all$avgb[1], color = "gray70") + 
    labs(y = "ATT estimate",
    x = "FY WIC EBT First Implemented",
    caption = paste0("ATT estimate across all groups is ",
                     as.character(round(all$avgb[1], digits = 3)),
                     " with a 95% confidence interval of (",
                     as.character(round(all$avgll[1], digits = 3)),
                     ", ",
                     as.character(round(all$avgul[1], digits = 3)),
                     ").")) +
    theme_light() + 
    scale_x_continuous(breaks = seq(2006, 2016, 2)) +
    theme(plot.caption = element_text(hjust = 0))
  
ggsave("csdid_tip_4zips_auth_all_group.png", plot = group_plot, width = 6.5, height = 4, units = "in")

group_plot <- ggplot(all, aes(x = gr_yr, y = b)) +
    geom_point(size = 3) +
    geom_ribbon(aes(ymin = all$avgll[1], ymax = all$avgul[1]),
              alpha = 0.25,
              color = "gray90") +
    geom_linerange(aes(ymin = ll, ymax = ul),
                alpha = 0.7,
                size = 0.8) +
    geom_smooth(method = lm, se = FALSE) +
    stat_poly_eq(formula = y ~ x,
                output.type = "numeric",
                mapping = aes(label = sprintf("Slope: %s",
                            c(formatC(round(stat(coef.ls)[[1]][[2, "Estimate"]], digits = 3))))),
                label.x = 3.5,
                label.y = -0.05) +
    geom_hline(yintercept = 0, lty = 2) +
    geom_hline(yintercept = all$avgb[1], color = "gray70") + 
    labs(y = "ATT estimate",
    x = "FY WIC EBT First Implemented",
    caption = paste0("ATT estimate across all groups is ",
                     as.character(round(all$avgb[1], digits = 3)),
                     " with a 95% confidence interval of (",
                     as.character(round(all$avgll[1], digits = 3)),
                     ", ",
                     as.character(round(all$avgul[1], digits = 3)),
                     ").")) +
    theme_light() + 
    scale_x_continuous(breaks = seq(2006, 2016, 2)) +
    theme(plot.caption = element_text(hjust = 0))
  
ggsave("csdid_tip_4zips_auth_all_group_fitline.png", plot = group_plot, width = 6.5, height = 4, units = "in")

# Define the colors for each chain status
colors <- brewer.pal(3, "Set2")[1:2]

# Define the shapes for each chain status
shapes <- c(16, 17)

chainindep <- filter(picture, controls == "no" & (subsample == "Chain" | subsample == "Independent")) %>%
    group_by(subsample) 
# Add new columns for color and shape mapping
chainindep$color <- factor(chainindep$subsample)
chainindep$shape <- factor(chainindep$subsample)

# Create the chain/independent plot
chainindep_plot <- ggplot(data = chainindep, aes(x = gr_yr, y = b, color = color, shape = shape)) +
    geom_linerange(aes(ymin = ll, ymax = ul),
                   alpha = 0.7,
                   size = 0.8) +
    geom_point(size = 3) +
    geom_hline(yintercept = 0, lty = 2) +
    scale_color_manual(values = colors, name = "Chain status") +
    scale_shape_manual(values = shapes, name = "Chain status") +
    geom_line(data = chainindep, aes(y = avgb, x = gr_yr, color = color)) + 
    geom_ribbon(data = chainindep, aes(ymin = avgll, ymax = avgul, color = color),
                   alpha = 0.1) +
    labs(y = "ATT estimate", x = "FY WIC EBT First Implemented", 
    caption = paste0("ATT estimate for the chain subsample is ",
                as.character(round(filter(chainindep, subsample == "Chain")$avgb[1] , digits = 3)),
                " with a 95% confidence interval of (",
                as.character(round(filter(chainindep, subsample == "Chain")$avgll[1], digits = 3)),
                ", ",
                as.character(round(filter(chainindep, subsample == "Chain")$avgul[1], digits = 3)),
                ").\nATT estimate for the independent subsample is ",
                as.character(round(filter(chainindep, subsample == "Independent")$avgb[1], digits = 3)),
                " with a 95% confidence interval of\n(",
                as.character(round(filter(chainindep, subsample == "Independent")$avgll[1], digits = 3)),
                ", ",
                as.character(round(filter(chainindep, subsample == "Independent")$avgul[1], digits = 3)),
                ").")) +
    theme_light() +
    theme(plot.caption = element_text(hjust = 0), legend.position="none") + 
    scale_x_continuous(breaks = seq(2006, 2016, 2)) +
    facet_wrap(~ subsample, ncol = 2)
  
ggsave("csdid_tip_4zips_auth_chainindep_group.png", plot = chainindep_plot, width = 6.5, height = 4, units = "in")


# Create the chain/independent plot with fitline
chainindep_plot <- ggplot(data = chainindep, aes(x = gr_yr, y = b, color = color, shape = shape)) +
    geom_linerange(aes(ymin = ll, ymax = ul),
                   alpha = 0.7,
                   size = 0.8) +
    geom_point(size = 3) +
    geom_smooth(method = lm, se = FALSE) +
    stat_poly_eq(formula = y ~ x,
                output.type = "numeric",
                mapping = aes(label = sprintf("Slope: %s",
                            c(formatC(round(stat(coef.ls)[[1]][[2, "Estimate"]], digits = 3)),
                              formatC(round(stat(coef.ls)[[2]][[2, "Estimate"]], digits = 3))))),
                label.x = 3.5,
                label.y = -0.05) +
    geom_hline(yintercept = 0, lty = 2) +
    scale_color_manual(values = colors, name = "Chain status") +
    scale_shape_manual(values = shapes, name = "Chain status") +
    geom_line(data = chainindep, aes(y = avgb, x = gr_yr, color = color)) + 
    geom_ribbon(data = chainindep, aes(ymin = avgll, ymax = avgul, color = color),
                   alpha = 0.1) +
    labs(y = "ATT estimate", x = "FY WIC EBT First Implemented", 
    caption = paste0("ATT estimate for the chain subsample is ",
                as.character(round(filter(chainindep, subsample == "Chain")$avgb[1] , digits = 3)),
                " with a 95% confidence interval of (",
                as.character(round(filter(chainindep, subsample == "Chain")$avgll[1], digits = 3)),
                ", ",
                as.character(round(filter(chainindep, subsample == "Chain")$avgul[1], digits = 3)),
                ").\nATT estimate for the independent subsample is ",
                as.character(round(filter(chainindep, subsample == "Independent")$avgb[1], digits = 3)),
                " with a 95% confidence interval of\n(",
                as.character(round(filter(chainindep, subsample == "Independent")$avgll[1], digits = 3)),
                ", ",
                as.character(round(filter(chainindep, subsample == "Independent")$avgul[1], digits = 3)),
                ").")) +
    theme_light() +
    theme(plot.caption = element_text(hjust = 0), legend.position="none") + 
    scale_x_continuous(breaks = seq(2006, 2016, 2)) +
    facet_wrap(~ subsample, ncol = 2)
  
ggsave("csdid_tip_4zips_auth_chainindep_group_fitline.png", plot = chainindep_plot, width = 6.5, height = 4, units = "in")